{
 "cells": [
  {
   "cell_type": "code",
   "id": "initial_id",
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "end_time": "2025-08-17T10:51:54.099861Z",
     "start_time": "2025-08-17T10:51:54.095961Z"
    }
   },
   "source": [
    "from scipy.integrate import odeint\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题\n",
    "plt.rc('font', size=10)\n",
    "plt.rc('font', family='SimHei')"
   ],
   "execution_count": 58,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T10:51:54.612185Z",
     "start_time": "2025-08-17T10:51:54.607885Z"
    }
   },
   "cell_type": "code",
   "source": [
    "m_A = 4866\n",
    "m_B = 2433\n",
    "f = 1760\n",
    "delta_m = 1091.099\n",
    "omiga = 1.9802\n",
    "k_3 = 528.5018  # 兴波\n",
    "k_5 = 80000  # 弹簧劲度系数\n",
    "k_6 = 10000  # 阻尼系数\n",
    "ro = 1025\n",
    "g = 9.8\n",
    "s = np.pi\n",
    "l_0 = 0.5  # 弹簧原长\n",
    "I_1 = 8289.43  # 浮子转动惯量\n",
    "h_2 = 0.5  # 振子的高度\n",
    "delta_I = 7001.914\n",
    "k_7 = 1655.909  # 兴波\n",
    "k_8 = 250000  # 弹簧劲度系数\n",
    "k_9 = 1000  # 阻尼系数\n",
    "k_10 = 8890.7  # 静水恢复力矩\n",
    "L = 2140\n",
    "\n",
    "# m_A = 4866\n",
    "# m_B = 2433\n",
    "# f = 4890\n",
    "# delta_m = 1165.992\n",
    "# omiga = 2.2143\n",
    "# k_3 = 167.8395  # 兴波\n",
    "# k_5 = 80000  # 弹簧劲度系数\n",
    "# k_6 = 10000  # 阻尼系数\n",
    "# ro = 1025\n",
    "# g = 9.8\n",
    "# s = np.pi"
   ],
   "id": "46ec89d2bacd559c",
   "execution_count": 59,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T10:51:55.268727Z",
     "start_time": "2025-08-17T10:51:55.265277Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 在SA函数外部初始化记录容器\n",
    "history = {\n",
    "    'temperature': [],\n",
    "    'best_P': [],\n",
    "    'current_P': [],\n",
    "    'k6': [],\n",
    "    'k9': [],\n",
    "}"
   ],
   "id": "fa504dcd8171dd89",
   "execution_count": 60,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T10:53:08.755543Z",
     "start_time": "2025-08-17T10:53:08.750105Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def plot_SA_history():\n",
    "    plt.figure(figsize=(15, 8))\n",
    "\n",
    "    # 子图1：损失函数变化\n",
    "    plt.subplot(2, 2, 1)\n",
    "    plt.plot(history['best_P'], 'r-', label='Best P')\n",
    "    plt.plot(history['current_P'], 'b--', alpha=0.5, label='Current P')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('P')\n",
    "    plt.title('P变化曲线')\n",
    "    plt.legend()\n",
    "    plt.grid(True, which=\"both\", ls=\"--\")\n",
    "\n",
    "    # 子图2：温度衰减曲线\n",
    "    plt.subplot(2, 2, 2)\n",
    "    plt.plot(history['temperature'], 'g-')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('Temperature')\n",
    "    plt.title('温度下降曲线')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    # 子图3：k6参数演化\n",
    "    plt.subplot(2, 2, 3)\n",
    "    plt.plot(history['k6'])\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('k6参数值')\n",
    "    plt.title('k6参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    # 子图3：k9参数演化\n",
    "    plt.subplot(2, 2, 4)\n",
    "    plt.plot(history['k9'])\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('k9参数值')\n",
    "    plt.title('k9参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "ee66a50ae1eeb76d",
   "execution_count": 64,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T10:53:10.157696Z",
     "start_time": "2025-08-17T10:53:09.827923Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def P(k_6, k_9):\n",
    "\n",
    "            # 定义一个方程组（微分方程组）\n",
    "        # 计算振子的转动惯量\n",
    "    def I_2(x_1, x_2):\n",
    "        delta_x = m_B * g / k_5\n",
    "        l_1 = x_2 - x_1 + l_0 - delta_x\n",
    "        l_2 = l_1 + h_2\n",
    "        return 1 / 3 * m_B / h_2 * (l_2 ** 3 - l_1 ** 3)\n",
    "\n",
    "\n",
    "    # 定义一个方程组（微分方程组）\n",
    "    def pfun(y, t):\n",
    "        y0, y1, y2, y3, y4, y5, y6, y7 = y\n",
    "        return np.array([\n",
    "            1 / (m_A + delta_m) * (-k_3 * y0 - ro * g * s * y1 + f * np.cos(omiga * t) + k_5 * (y3 - y1) + k_6 * (y2 - y0)),\n",
    "            y0,\n",
    "            1 / m_B * (-k_5 * (y3 - y1) - k_6 * (y2 - y0)),\n",
    "            y2,\n",
    "            1 / (I_1 + delta_I) * (-k_7 * y4 - k_10 * y5 + L * np.cos(omiga * t) + k_8 * (y7 - y5) + k_9 * (y6 - y4)),\n",
    "            y4,\n",
    "            1 / I_2(y1, y3) * (-k_8 * (y7 - y5) - k_9 * (y6 - y4)),\n",
    "            y6,\n",
    "        ])\n",
    "\n",
    "\n",
    "    t = np.arange(0, 1000, 0.2)  # 创建自变量序列\n",
    "    soli = odeint(pfun, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], t)  # 求数值解\n",
    "    # print(soli[:, 0])\n",
    "    # plt.figure(figsize=(5, 2))\n",
    "    # plt.plot(t, soli[:, 1], 'b', label=\"浮子位移\")\n",
    "    # plt.plot(t, soli[:, 3], 'r', label=\"振子位移\")\n",
    "    # plt.legend()\n",
    "    # plt.show()\n",
    "    # print(soli[:, 6])\n",
    "\n",
    "    re = 0\n",
    "    for i in range(200, len(soli[:, 2])):\n",
    "        re += (k_6 * (soli[:, 0][i] - soli[:, 2][i])**2 + k_9 * (soli[:, 6][i] - soli[:, 4][i])**2)*0.2\n",
    "\n",
    "\n",
    "    return re/960\n",
    "re = P(40000, 98227)\n",
    "print(re)"
   ],
   "id": "a93479dfcaa0cc8b",
   "execution_count": 65,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T10:54:01.974597Z",
     "start_time": "2025-08-17T10:53:11.488456Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def SA(iter, t0, tf, alpha):\n",
    "    t = t0\n",
    "    k6c = 20000\n",
    "    k9c = 20000\n",
    "    Pc = P(k6c, k9c)\n",
    "    k6b = k6c\n",
    "    k9b = k9c\n",
    "    Pb = Pc\n",
    "    for i in range(iter):\n",
    "        k6n = k6c + np.random.randint(-50000, 50000)\n",
    "        k6n = np.clip(k6n, 0, 100000)\n",
    "        k9n = k9c + np.random.randint(-50000, 50000)\n",
    "        k9n = np.clip(k9n, 0, 100000)\n",
    "        Pn = P(k6n, k9n)\n",
    "        if Pn > Pc or np.random.rand() < np.exp((Pn-Pc)/t):\n",
    "            k6c = k6n\n",
    "            k9c = k9n\n",
    "            Pc = Pn\n",
    "            if Pc > Pb:\n",
    "                k6b = k6c\n",
    "                k9b = k9c\n",
    "                Pb = Pc\n",
    "        t = t*alpha\n",
    "        if t < tf:\n",
    "            break\n",
    "        print(f\"第{i}轮,kb:{k6b},Pb:{Pb}, kc:{k6c},Pc:{Pc}\")\n",
    "        history['temperature'].append(t)\n",
    "        history['best_P'].append(Pb)\n",
    "        history['current_P'].append(Pc)\n",
    "        history['k6'].append(k6c)\n",
    "        history['k9'].append(k9c)\n",
    "    return k6b, k9b, Pb\n",
    "re = SA(800, 2000, 0.001, 0.9)\n",
    "print(re)\n",
    "plot_SA_history()"
   ],
   "id": "5230104b2da26a12",
   "execution_count": 66,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T10:50:17.315794Z",
     "start_time": "2025-08-17T10:50:17.311706Z"
    }
   },
   "cell_type": "code",
   "source": [
    "for key in history.keys():\n",
    "    history[key].clear()\n",
    "print(history['temperature'])"
   ],
   "id": "c1670b45ca025376",
   "execution_count": 57,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T08:55:47.470371Z",
     "start_time": "2025-08-17T08:55:47.467594Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 在SA函数外部初始化记录容器\n",
    "history2 = {\n",
    "    'temperature': [],\n",
    "    'best_P': [],\n",
    "    'current_P': [],\n",
    "    'k': [],\n",
    "    'alpha': []\n",
    "}"
   ],
   "id": "99ec20b8d5701ec2",
   "execution_count": 13,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T08:56:06.336690Z",
     "start_time": "2025-08-17T08:56:06.331558Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def plot_SA_history2():\n",
    "    plt.figure(figsize=(10, 8))\n",
    "\n",
    "    # 子图1：损失函数变化\n",
    "    plt.subplot(2, 2, 1)\n",
    "    plt.plot(history2['best_P'], 'r-', label='Best P')\n",
    "    plt.plot(history2['current_P'], 'b--', alpha=0.5, label='Current P')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('P')\n",
    "    plt.title('P变化曲线')\n",
    "    plt.legend()\n",
    "    plt.grid(True, which=\"both\", ls=\"--\")\n",
    "\n",
    "    # 子图2：温度衰减曲线\n",
    "    plt.subplot(2, 2, 2)\n",
    "    plt.plot(history2['temperature'], 'g-')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('Temperature')\n",
    "    plt.title('温度下降曲线')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    # 子图3：k参数演化\n",
    "    plt.subplot(2, 2, 3)\n",
    "    plt.plot(history2['k'])\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('k参数值')\n",
    "    plt.title('k参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1))\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "        # 子图4：alpha参数演化\n",
    "    plt.subplot(2, 2, 4)\n",
    "    plt.plot(history2['alpha'], c='orange')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('alpha参数值')\n",
    "    plt.title('alpha参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1))\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "398b619d4091fcaf",
   "execution_count": 15,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T08:56:08.899717Z",
     "start_time": "2025-08-17T08:56:08.786751Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def P2(k, alpha):\n",
    "        # 定义一个方程组（微分方程组）\n",
    "    def pfun(y, t):\n",
    "        y0, y1, y2, y3 = y\n",
    "        return np.array([\n",
    "            1 / (m_A + delta_m) * (-k_3 * y0 - ro * g * s * y1 + f * np.cos(omiga * t) + k_5 * (y3 - y1) + k*(abs(y2-y0))**alpha * (y2 - y0)),\n",
    "            y0,\n",
    "            1 / m_B * (-k_5 * (y3 - y1) - k*(abs(y2-y0))**alpha * (y2 - y0) ),\n",
    "            y2\n",
    "        ])\n",
    "    t = np.arange(0, 500, 0.2)  # 创建自变量序列\n",
    "    soli = odeint(pfun, [0, 0, 0, 0], t)\n",
    "    # print(soli[:, 0])\n",
    "    # plt.figure(figsize=(5, 2))\n",
    "    # plt.plot(t, soli[:, 1], 'b', label=\"浮子位移\")\n",
    "    # plt.plot(t, soli[:, 3], 'r', label=\"振子位移\")\n",
    "    # plt.legend()\n",
    "    # plt.show()\n",
    "\n",
    "    re = 0\n",
    "    for i in range(200, len(soli[:, 2])):\n",
    "        re += k*(abs(soli[:, 0][i] - soli[:, 2][i]))**alpha * (soli[:, 0][i] - soli[:, 2][i])**2*0.2\n",
    "    return re/460\n",
    "re = P(40000)\n",
    "print(re)"
   ],
   "id": "a6ae873eb41a9ad",
   "execution_count": 17,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T08:57:09.475421Z",
     "start_time": "2025-08-17T08:56:09.609573Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def SA(iter, t0, tf, alpha):\n",
    "    t = t0\n",
    "    kc = 50000\n",
    "    ac = 0.5\n",
    "    Pc = P2(kc, ac)\n",
    "    kb = kc\n",
    "    ab = ac\n",
    "    Pb = Pc\n",
    "    for i in range(iter):\n",
    "        kn = kc + np.random.randint(-50000, 50000)\n",
    "        kn = np.clip(kn, 0, 100000)\n",
    "        an = ac + np.random.normal(0,0.5)\n",
    "        an = np.clip(an, 0, 1)\n",
    "        Pn = P2(kn, an)\n",
    "        if Pn > Pc or np.random.rand() < np.exp((Pn-Pc)/t):\n",
    "            kc = kn\n",
    "            Pc = Pn\n",
    "            ac = an\n",
    "            if Pc > Pb:\n",
    "                ab = ac\n",
    "                kb = kc\n",
    "                Pb = Pc\n",
    "        t = t*alpha\n",
    "        if t < tf:\n",
    "            break\n",
    "        print(f\"第{i}轮,kb:{kb},ab:{ab}, Pb:{Pb}, kc:{kc},ac:{ac},Pc:{Pc}\")\n",
    "        history2['temperature'].append(t)\n",
    "        history2['best_P'].append(Pb)\n",
    "        history2['current_P'].append(Pc)\n",
    "        history2['k'].append(kc)\n",
    "        history2['alpha'].append(ac)\n",
    "    return kb, ab, Pb\n",
    "re = SA(1000, 2000, 0.001, 0.95)\n",
    "print(re)\n",
    "plot_SA_history2()"
   ],
   "id": "273ca2ca10762f6b",
   "execution_count": 18,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-17T04:00:44.598253Z",
     "start_time": "2025-08-17T04:00:44.594343Z"
    }
   },
   "cell_type": "code",
   "source": [
    "for key in history2.keys():\n",
    "    history2[key].clear()\n",
    "print(history2['temperature'])"
   ],
   "id": "e707762b7dd7ce94",
   "execution_count": 141,
   "outputs": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
